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Abstract 

The liquid-gas phase diagram for polydisperse dipolar hard-sphere fluid with polydisper- 
sity in the hard-sphere size and dipolar moment is calculated using extension of the recently 
proposed thermodynamic perturbation theory for central force (TPT-CF) associating po- 
tential. To establish connection with the phase behavior of ferro colloidal dispersions it is 
assumed that the dipole moment is proportional to the cube of the hard-sphere diameter. 
We present the full phase diagram, which includes cloud and shadow curves, binodals and 
distribution functions of the coexisting daughter phases at different degrees of the system 
polydispersity. 
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I. INTRODUCTION 



In this paper we consider the liquid-gas phase behavior of polydisperse dipolar 
hard-sphere mixture. Recently, liquid-gas phase equilibria in monodisperse dipo- 
lar hard-sphere fluid, Yukawa dipolar hard-sphere fluid and Shtockmayer fluid were 
studied using thermodynamic perturbation theory for central force (TPT-CF) as- 
sociating potential . In this study we propose extension of the TPT-CF, 
which enables us to investigate the phase behavior of polydisperse mixture of the 
dipolar hard spheres with polydispersity in both hard-sphere size and dipole mo- 
ment. We call this extension as extended TPT-CF (ETPT-CF). In a certain sense 
ETPT-CF can be seen as a combination of Wertheim's TPT 5|, |6| for associating 
fluid with association taking place due to off-center attractive sites and TPT-CF, 
which allows multiple bonding of each site. In our theory we have several Wertheim's 
type of associating sites with the possibility for each site to be bonded to several 
other sites (in Wertheim's TPT each site is only singly bondable). Final expressions 
for thermodynamical properties of polydisperse dipolar hard-sphere fluid is written 
in terms of the finite number of distribution function moments, i.e. in the frame- 
work of ETPT-CF this system belongs to the family of the so-called truncatable 
free energy models. This property enables us to calculate the full liquid-gas phase 
diagram (including cloud and shadow curves and binodals) and to study effects of 
fractionation on the level of the distribution functions of coexisting daughter phases. 



II. EXTENDED THERMODYNAMIC PERTURBATION THEORY FOR 
CENTRAL FORCE ASSOCIATIVE POTENTIAL 

A. Analysis and classification of diagrams 

We consider a multicomponent fluid mixture consisting of n species with a number 
density p = Pa & t a temperature T (/3 = 1/ksT), where p a is the density of the 
particles of a species. The particles of the species a and h interact via the pair 
potential [7^(12), which can be written as a sum of the reference U^(12) and 
associative U^ s (12) parts 
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U ab {12) = U£ f {12) + U£,{12), (1) 

where 1 and 2 denote positions and orientations of the particles 1 and 2. We 
assume that associative part of the potential can be represented as a sum of M a x M& 
terms, i.e. 

E# a (12) = 5>^(12), (2) 

where the lower indices K and L take the values A, B,C, . . . and A, B,C, . . ., 

respectively. These values specify splitting of the total associating potential U^ s (12) 
into several particular pieces. For example in the case of the models utilized by 
Wertheim 6| these indices denote off-center attractive sites and in the case of the 
Mersedes-Benz (MB) type of the models \7] or cone models [8| they stand for the 
type of the hydrogen bonding arms. Hereafter we will refer to these indices as to the 
site indices, keeping in mind that they may have more general meaning. Here M a 
and Mf, are the number of such sites on the particles of a and b species, respectively. 
According to (pQ) and ([2} Mayer function f ab (12) for the total potential (p]) takes the 
following form: 



U(12) = f r %(12) + ef ef (12) \H[1 + /&(12)] - 1 \ , (3) 
where we are using the usual notation: 



. kl 



e(12) = exp [-/3C/(12)], /(12) = e(12) - 1, (4) 

For the sake of diagrammatic analysis we will follow Wertheim 6| and instead of 
circles we introduce hypercircles to represent particles in diagrammatic expansions. 
Each hypercircle is depicted as a large open circle with small circles inside denoting 
the sites. Corresponding cluster integrals are represented by the diagrams build on 
a hypercircles connected by f re f and e re f bonds and a site circles connected by the 
associating bonds Jkl- Due to the decomposition of the Mayer function / a f,(12) (j3J) 
we will have the following diagrammatic expressions for the logarithm of the grand 
partition function H and for the one-point density p a (l) in terms of the activity z: 
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InS = sum of all topologically distinct connected diagrams consisting of field 
z hypercircles , f re f, e re / and fxL bonds. Each bonded pair of z hypercircles has 
either f re f, or e re f and one or more fxi bonds. 

p a (l) = sum of all topologically distinct connected diagrams obtained from 
InS by replacing in all possible ways one field z hypercircle by a z a (l) circle labeled 
1. 

Here z a {i) = z a exp [— f3U a (i)\, i denotes position and orientation of the parti- 
cle i, and U a (i) is an external field. For a uniform system z n (l) = z„. Following 

nnnn 

[11, 12|, [5|, [6|] we introduce the definition of the s-mer diagrams. These are the 
diagrams consisting of s hypercircles, which all are connected by the network of fxL 
bonds. The site circles, which are incident with more than m a K f^ L bonds are called 
oversaturated site circles. We consider now the set of oversaturated site circles with 
each pair connected by at least one path formed by the circles from the same set. 
The subdiagram involving this set of circles, together with the site circles adjacent 
to them and fxL bonds connecting all these circles , we call the oversaturated 
subdiagram. The set of all possible s-mer diagrams can be constructed in three 
steps: (i) generating the subset of all possible connected diagrams with only fxL 
bonds, (ii) inserting combined bond e re f = f re f + 1 between all pairs of hypercircles 
with the site circles, which belong to the same maximal oversaturated subdiagram 
and (iii) taking all ways of inserting a f re f bond between the pairs of hypercircles, 
which were not connected during the previous two steps. As a result the diagrams, 
which appear in InS and p(l), can be expressed in terms of the s-mer diagrams: 

InS = sum of all topologically distinct connected diagrams consisting of s- 
mer diagrams with s = l,...,oo and f re f bonds between pairs of hypercircles in 
distinct s-mer diagrams. 

The procedure for obtaining the expression for p a (l) from InS remains un- 
changed. 

The diagrams appearing in the z expansion of the singlet density p a (l) can be 
classified with respect to the number of f^ L bonds associated with the labeled z a {\) 
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hypercircle. We denote the sum of the diagrams with < m a K associating bonds 
connected to the site K (K — A, B,C, . . .), which belongs to the particle of species 
a as Bi B ,c io ,...(l)- Any s ^ e which is connected to %k > Tn K associating bons, 



will be denoted as K m a In what follows we will use also a condensed version of the 



notation, i.e. 



P% a ,b ib ,c^)=pU b ,cJ 1 ) = pU 1 ) (5) 
where {i} = ia^b, ic, ■ ■ ■■ The set {i} with all indices, except one index i K , equal 0, 
will be denoted as ix, i.e. {i} = 0, . . . , 0, ix, 0, . . . , = i^, so that for any quantity 
x^y we have 

X {i} = X 0,...,0,i K ,0,...,0 = x k 1r = x i K (6) 
Thus PaiX) can be written as follows 

m%m B ,... {m a } 

B. Topological reduction 

The process of switching from the activity to a density expansion goes in the 



same fashion as in Refs. 



11,33. However, to proceed it is convenient to use an 
operator form of notation. The operators are introduced in a manner similar to 
that presented in references to which we refer the reader for more details. 

We associate with each labeled I hypercircle an operator em.(i) with the following 
properties: 

e {i}(0 = 0, if any i K > m a K , 

e£ } (0 = 1, if all % K = 0, (8) 

where {i + j} = %a + Ja, + js, ic + jc, ■ ■ •• The one-point quantities, which, for 
convenience, are denoted by can be presented as illustrated below: 

{m a } 

x a (i) = E zh^Kyi 1 )- ( 9 ) 

{i}=0 
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The operators e|o are straightforward generalization of the operators introduced 
earlier nag. Thus, the rules of manipulation with the new quantities Xa are 
similar to that discussed before. In particular, the usual algebraic rules apply to 
these quantities and analytical functions of Xa £1X6 defined by the corresponding 
power series. Similar, as in references GIB , it is convenient to use the angular 
brackets to specify the operation 

(x a ) = X{ m a}- (10) 

In the case of several labeled circles the subscripts on the brackets denote the circle 
to which procedure ( ITOl) is to be applied. 

Analyzing the connectivity of the diagrams in p a (l) at a labeled z a (l) hypercircle 
we have 



p a (l)/z a (l) = exp[ca(l)], (11) 

where 0/^(1) with {i} ^ {0} denotes the sum of diagrams in p%y(l) / p'LJl) for which 
the labeled 1 hypercircle is not an articulation circle. Similarly C| |(l) denotes the 
sum of diagrams in P| j(l) / z a (l) for which hypercircle 1 is not an articulation circle. 
Elimination of the diagrams containing field articulation circles can be achieved by 
switching from an activity to a density expansion. To do so we adopt the following 
rule: each field hypercircle z a , with bonding state of its sites represented by the 
set {/}, in all irreducible diagrams c a is replaced by a a%y hypercircle, where jx = 
m °k — Ik (K — A,B,. . .) for m a K — l K > and j K = for m a K — l K < 0. The new 
quantities a^(l) are connected to the densities via the following relation: 

{m a } 

* (i) = p.(i) Y, e U 1 )- ( 12 ) 

{i}=0 

1 



This relation can be inverted expanding ^{I}=o e {i}(l) 



into a power series, i.e. 



=*«(i) n [i-^wi ( i3 ) 



K=A 



Now the diagrammatic expansions for c?^ can be expressed in terms of irreducible 
diagrams. To present this result in compact and convenient form we introduce a 
sum of the diagrams defined as follows: 
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= sum of all topologically distinct irreducible diagrams consisting of s- 
mer diagrams with s = l,...,oo and f re f bonds between pairs of hypercircles in 
distinct s-mer diagrams. All hypercircles are field circles carrying the a-factor 
according to the rule formulated above 

Functional differentiation of with respect to o-^ ma _^ gives an expression 
for c a {l} : 



5c<® 



^ 6a? 



(14) 



{m a —i} 



C. Extended thermodynamic perturbation theory for central force associ- 
ating potential 



Now we are in a position to rewrite the regular one-density virial expansion for 
the pressure P in terms of the density parameters <r a (l). Following the scheme, 
proposed earlier j^, |j, 5, 6] we have expression for the pressure in operator form 



-,(0) 



and explicitly 



0pv = y.J 



{m a } 

^ a { V-, } (l)4 } (l) 

{i}=0 



d(l)+c 



(o) 



(15) 



(16) 



where V is the volume of the system. Similarly, as in [5|, [6J one can verify that 
these expressions satisfy the regular thermodynamic relation p a = /3dP/dfi a , where 
p a = J Pa(l) d(l) and /i a is the chemical potential. This can be achieved by taking 
a variation of (I15p (or ([175]) ) and combining (fTTj) . (1131) and ([Til) . The corresponding 
expression for Helmholtz free energy is 



Pa(l)ln 



An 



{m a } 



^(l)-c (0) , (17) 



where A a is the thermal de Broglie wavelength. This expression is derived using the 
regular thermodynamic expression for Helmholtz free energy A = Yl a N a ^a — PV 
together with relation 



PNalla = / Pa {l) 



In — —r c 



\ 4o}( 1 ' 



d(l), 



(18) 



which follows from ffTTl) . written for P{ }- Here iV a is the number of particles of 
species a in the system. 

Helmholtz free energy in excess to its reference system value A re f is obtained by 
subtracting corresponding expression for A re f from ffTTj) . i.e. 



a J Pa{l) {^0 



/3 {A -Aref) 
{m a } 



{m°-i}V / {*} 



DcUD 



d(l) - (c (0) - C) , (19) 



where is the corresponding sum of the diagrams for the reference system. Or- 
dering the virial expansion (TT9"j) with respect to the number of associating f KL bonds 
and neglecting the terms with more than one associating bond we have 



.(o) _ Jo) 

C ref 



lY, / flS / (12)<^„fl)/.„,(12)^(2)> 12 r/(l)r/ t 2) 



ab 



and 



c a (l) - cfo(l) = £ / < f (12)(/ a6 (12)a a (2)) 2 d(2), 



(20) 



(21) 



where g^^(12) is the reference system distribution function and 

/ o6 (12) = 5>^(l)/^(12)4(2). (22) 

Due to single bond approximation c?^ = for all values of the set {i}, except 
for {i} = and {i} = i K with i K — 1. This property together with (ITTj) yield the 
following relations: 



Px 1 (l)/P{o}( 1 ^ 



and 



n 



pfo}(!: 



(23) 



for z K G {«} (24) 



The set of relations (1201 . (I2TI) and (|23|) defined all the quantities needed to calculate 
the Helmholtz free energy of the system f[T9"j) . provided that the properties of the 
reference system are known. 
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Finally it is worth noting, that ETPT-CF theory developed here reduces to TPT1 
proposed by Wertheim [6(, if for all sites single-bonding condition m a K = 1 is as- 
sumed. In the other limiting case of only one site per particle ETPT-CF will coincide 
with TPT-CF developed earlier [l, 2, 4]. 



D. Extended TPT-CF for two sites with double-bonding condition 

The theory presented in the previous section is quite general and can be applied 
to a number of different situations. However in the present study we are interested 
in the version of the theory for the model with two sites both of which can be bonded 
twice. More specifically we are interested in the extension and application of the 
theory to study the phase behavior of polydisperse dipolar hard-sphere mixture. 

We assume that each of the particles in the system have two doubly bondable 
attractive sites, A and B, i.e. we have: M a = 2 and m a A = m a B . We also assume, 
that attractive interaction is acting only between the sites of the same sort. Using 
these suggestions, relations ( fTTj) and ( |T2l and taking into account that the system 
is uniform, the density parameters (Ta 2 b 2 = Pa <t a 1 b 2 = an d (T a 2 b 1 = ^Bi a can 
be expressed in terms of c a Ki 

1 „ 



Pa 



1 a A B 



%-a 
Ai 



(X B 1 



i + K) 



2 (T A B K A 1 



2 (J A B K B 1 



(25) 

(26) 
(27) 



where K takes the values A and B and n a Ki = 1 + c a Ki . These two equations give 



and 



2pa^Xi 



-1 



1 1 + 



2 ' 



(28) 



(29) 



In turn, using (12ip . for k\ we have 



K 



T, 1 



ab ~a 
KK°Kx i 



(30) 
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where 



Ikk = J fl&(12)/&(12) d(2). (31) 
Combining ( |20|) . ( 1291) and ( 1301) expression for Helmholtz free energy (fl~9]) can be 
written in terms of k^- parameters 



A — A, 



ref 



V 



Pa 



(32) 



which satisfy the following set of equations: 



K 



rab 
~2 1 KK 



(33) 
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Chemical potential A/i Q and pressure AP in excess to their reference system 
values can be obtained using standard thermodynamical relations: 



Pa Pre/ 



d[(A-A ref )/V] 



P - P ref = ^ Pa (Pa - Pref) 



A- A 



ref 



• (34) 



df)a > ' -r*j ^r-v- rre;; y 

Finally, the average size of the clusters, which appear in the system, can be 
characterized by the average length of the chain Lk formed by either A-bonded 
(K = A) or 5-bonded (K = 5)particles. This quantity is defined by the following 
expression 



Ea ( a K,end + <,mid + <) 



(35) 



where a a K end is the fraction of singly i^-bonded particles (fraction of the chain ends), 
a Kmid i s the fraction of doubly i^-bonded particles (fraction of the chain middles) 
and «q is the fraction of nonbonded particles. For these fractions we have 



P a A,end = &Ai ~ a A B 2 , P a A,mid = P ~ a A V P a = a A B - 



(36) 



Substituting these expressions into expression for La (135]) and using ([28]) . (|29|) and 
expression for cr A B2 , 

l r 2 i 

a A B 2 = -jPaOA Q B a + lj , (37) 

we get the final expression for La in terms of K a Kl : 

4-(l-«&J 



J K = 2_^Pa 



4- 






>1 




1 + fef 









E 



Pa" 



-1 



1 + fa) 1 + fa) 



(38) 
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where if K = A then K = B and if K = B then K = A. 
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III. LIQUID-GAS PHASE BEHAVIOR OF POLYDISPERSE DIPOLAR 
HARD-SPHERE FLUID 

A. The model 

We consider polydisperse dipolar hard-sphere fluid mixture with a number den- 
sity p and polydispersity in both hard-sphere diameter a and dipolar moment d^. 
We assume, that the dipole moment is proportional to the particle volume, i.e. 
dfj, ~ a 3 . Thus the type of the particle is completely defined by its hard-sphere 
size and hereafter we will be using a instead of the indices a,b, . . . to denote the 
particle species. We also assume, that hard-sphere size of the particles is distributed 
according to normalized distribution function F(a) > 0, 

oo 

F(a) da = 1. (39) 

Interaction between particles of species o\ and 2 in our system is described by the 
following pair potential: 

U{r,a x a 2 ) = U hs {r,a x a 2 ) + U M (12, a l a 2 ), (40) 

where Uh s (r, cr\cr 2 ) is the hard-sphere potential and UM l (r,<j\<j<i) is the dipole-dipole 
potential, given by 

dn((7\)dn{(72) 

(7^(12, a 1 a 2 ) = -3 [2 cos y?i cos y9 2 - sin ipt sin tp 2 cos ((pi - 2 )J • (41) 

Here (fx and <p 2 denote the angles between the dipole vectors and vector, which 
joins the centers of the two particles, and <pi and <p 2 are the azimuthal angles about 
this vector. To proceed we have to split the total potential ( l4T)i) into the refer- 
ence and associative pieces. We assume, that the reference part of the potential is 
represented by the hard-sphere part Uh s (r, <7i<7 2 ) and associative part by the dipole- 
dipole potential Udd(r, a \a 2 ). At the contact distance 0\ 2 = (01 + cr 2 )/2 the latter 
potential has two equal potential minima of the depth — 2^(01 )d M (02) / af 2 at "nose- 
to-tail" configuration (pi = ip 2 = 0, ipi = tp 2 = n. These minima are responsible 
for the formation of the chains of the particles in the system. In addition there 
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are twice less deep minima (— d^(ai)d^(a 2 )/af 2 ) at antiparallel configuration with 
(pi — <f 2 — 7r/2, </>i — 02 = 7T- The latter minima cause the formation of the 
network connecting the chains. According to earlier theoretical and computer simu- 



lation studies 
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111 ] competition between chain formation and network formation 



defines the existence of the liquid-gas phase transition in the dipolar hard-sphere 
fluid. To account for this effect we propose the following splitting of the total asso- 
ciative potential U ass (12, a x a 2 ) = £7^(12, cr\a 2 ): 



U BB (12,a 1 cT 2 ) = e(ip 1 )e(<p 2 ) U M (12, 0102), (42) 

U AA (12, 0x0 2 ) = [1 - fa) 6 fa)] U dd (12, 0!0 2 ), (43) 

where Qfa = H (ir/2 + <p — ip) H {ti/2 — ip + <p) and H(x) is the Heaviside step 
function. Here ipo plays a role of the potential splitting parameter. For ip = n/2 
we have that £7aB(12, 0i0 2 ) = £^(12, 0-102) and Uaa(^2, aia 2 ) = 0. On the other 
hand ip = gives: L r s(12,0 1 2 ) = and £74^(12,0102) = £^m(12, 0102)- In both 
limiting cases the theory developed will treat the system as a polydisperse mixture 
of the hard-sphere chains. For the intermediate values of ipo the energy minima 
at "nose-to-tail" configuration are included into Uaa{^-2, <Ji<j 2 ) and network forming 
minima appear in £7^5(12,0102). Our choice for the particular value of (po will be 
discussed later. 



B. Thermodynamic properties 

For a general multicomponent dipolar hard-sphere mixture thermodynamic prop- 
erties can be obtained using solution of the set of nonlinear equations ( 1331) and 
expression for Helmholtz free energy fl32|) . However even for multicomponent case 
solution of this equation rapidly becomes involved as the number of components 
increases. As we proceed to polydisperse case solution of the polydisperse version 
of equation (J33l becomes intractable, since now we have to deal with the following 
integral equation: 

k k {(Ji) = 1 + 2p / F(a 2 ) — — 2—- da 2 , (44) 

JO l + K K\ a V 

where we have dropped the lower index 1, i.e. /^(o) = Kk{o). In order to solve 
this equation we propose here to interpolate the key quantity of the theory, the 
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volume integral I kl{^ 1^2)1 using a sum of Ny Yukawa terms. Since the reference 
system pair distribution function g re f(r,aia 2 ) is independent of mutual orientation 
of the particles for the integral (13T|) we have 



POO 

Ikk{viV2) =471 r 2 g re f(r, <J\ct 2 ) fia((r, a x a 2 ) dr, 
Jo 



(45) 



where fKK{f,(Jia 2 ) is orientation averaged Mayer function for associative potential 
C/k-jc(12, era). We assume, that ]kk{j'- i V\V'i) can be represented in the following 
form: 

fKK{r,o,o 2 ) = ±£A%\* 1 )A%\* a )e-®^\ (46) 

n 

Parameters Ajp(a) and ZfP are obtained using the interpolation scheme, which is 
presented and discussed in the Appendix A. Using (1431) and (j4"6"|) . we have 

W<7i<7 2 ) = ^ ^(^^(^G^V^), (47) 

n 

where G^\a\a 2 ) = re _z ^ >r 'g re f(r, <Jia 2 )dr is the Laplace transform of the radial 
distribution function g re f(r, aia 2 ). We will be using here Percus-Yewick approxima- 
tion for hard-sphere radial distribution function, since analytical expression for its 
Laplace transform is known 



,0) 



A 



y (n)\2 nW 
-A" J K 



zy [a 12 + crta 2 —m 2 J + 1 + 7^3+ 



where 



(n) 

, ( («) o 00 , («) 



n (n) A 2 2?r / i , 1 A ( (») , X 

= A _ 7h ( 1 + 2 7rms J 1 m *>o + 2 m2 



(4£ 



2« ^Am&i + ^ir 



m S ( m 2 + 2ms,) - ( m j?i 



(49) 



(m 2 + 2m\ 

A = 1 - 7rm 3 /6. (50) 

Here m; are the moments and are the generalized moments of the distribution 
function F(cr). Symbolically expression for these moments can be represented as 
follows: 

POO 

m = m{a)F{a) da. (51) 
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Hereafter all the quantities denoted as m with certain set of indices will represent 
the generalized moments defined by IpTD . Corresponding expressions for m(a) are 
collected in the Appendix B. Inserting (|47|) into equation (jUJ) and using (115]) . we 
find 



TVy 2 

n j=l 



(52) 



where C^- satisfies the following set of equations: 



C 



(n) 

AJ 



2 / ^F((r)^((r)- 



r(0 



,(0 



1 2 



(53) 



Here 



_ J[_ (n) 



p. 



A', 2 



2z v; D 

ZZ K K 
71 

o~( n ) n( n ) 

ZZ K K 



1 



(n) 

-m 2 + 77% o 



P 



(n) 



A 



if,3 



» 



(n) 
K 



2A V ms ^ m ^2 



(54) 
(55) 

(56) 

(57) 



Thus solution of the integral equation ( j44l) for the unknown function kk{&) now 
is reduced to solution of a set of equations for 4iVy unknown constants C^. This 
solution can be used to calculate % (cr), which in turn can utilized to calculate ther- 
modynamical properties of the system via Helmholtz free energy (l32l . Generalizing 
expression for Helmholtz free energy (l32l) for polydisperse system, we have 

B B 



A- A 



ref 



V 



fOO 

n I F(a) 

'0 



KAcr 



i)-£ 



K K {a\ 



K=A 



K=A 



1 - k k (g) 
l + K 2 K {a) 



da. 



(58) 



Now we can use the standard relation between Helmholtz free energy and chemical 
potential (1M|) . generalized to polydisperse case 



(59) 



where #/(!> {P(<r)} denote functional differentiation with respect to the distribution 
F(a). We find 

R R 



1 B 73 

(3 \p(a) - nref{a)\ = ^ - In - ]J + l) - £ 
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3 x pH 2 



v 



3=1 



(n) 0L/ K,j 

K ' j F~(a) 



Here 



(n) 

/%,2 



I ( m (»>°) _ n r (n A r {n) 



(n) 
Pk,3 



K.2 



m 



(n,0) 

jo 



P U K,2 



- m 
2 

(n,0) 
Jf,0 



(n,0) 



(n)' 



(n) \ p (n) 



>,0) 



(n) \ p (n) 



(60) 

(61) 
(62) 
(63) 
(64) 



functional derivatives SP^j/SF(a) and SD^' /5F(a) are presented in the Appendix 



(n) 



B and functional derivatives 5C^/8F(a) can be obtained from the solution of the 
set of linear equations, which follows from (1531) upon its functional differentiation 
with respect to F(cr), i.e. 



2 SC (K) 



1,2, K 



A,B, 



(65) 



where 
.(JO 



5Cj /8F(a) = 5C\"'J5F(a) and the elements of the matrices and 



t(n) 



(*0 



Rj- (cr) are collected in the Appendix B. 

The pressure expression follows from (l3"4"|) . generalized to polydisperse case 



P - P 



re/ — P 



F(a) [pi(a) - fi re f(cr)] da - 



A- A 



ref 



V 



(66) 



Using this expression together with expression for the chemical potential fl60l) . we 
find 



m, 



N Y B / 3 2 
n if+A \j"=l j=l 



(67) 



where 



J 



(n) 



v ' da, 



F(a) 



5F(a) 

r(n) 



o(n) 



F(a) 



(n) 
K,j 



5F(a) 



da. 



(6* 



(n) 



Expression for the integral J^-j- is presented in the Appendix B and integral S K 
can be obtain from the solution of the set of linear equations 



£M<«>Sf> 

3=1 



E 



(JO 



1,2, K 



A,B, 



(69) 



which follows from the set of equations ( 1651) . Here 
of the matrix are presented in the Appendix B. 



Sftj and the elements 
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Expressions for the chemical potential f[6"0l and pressure (I6T1) are the final ex- 
pressions to be used in the phase equilibrium calculations. The properties of the 
reference system (chemical potential /i re /(cr) and pressure P re f) are described here 
using polydisperse versions of the Mansoori et al. [13j expressions 

2 



,(ez), 



(J 



7ra [ 1 



m 3 J 



+ 



2A 

(3P re f 



—a 

3 



1 + A 



m 3 



7T 



In A + m 2 1 + 



77I3A 



m 



A 



A 



+ T m 2 7T m l + o 



1 ml 



1 

A 



7T 



P + ^^1^2 



7T 



12A 



2 m 2- 



3 7B3A 

7T\3 1 



6/ A 



2 m 3 m 2 



(70) 



(71) 



where p^/ ( a ) * s ^^e reference system chemical potential in excess to its ideal gas 



value. 



C. Phase equilibrium conditions 

One can easily see that the thermodynamical properties of the model at hand 
obtained above are defined by the set of the finite number of the distribution function 
moments, i.e., four regular moments (mi, I = 0, 1, 2, 3) and l + 10iVy + 3iVy ( Ny + 1) 
generalized moments (m M , m^\ m^™' '; i = 0, 1,2; j = 0, 1; / = 0, 1,2; if = 
A, B). Note, that m^f 1 = m^ 1 ^. Thus the polydisperse mixture of dipolar hard 
spheres treated within ETPT-CF belongs to the class of truncatable free energy 
models jl4j]. This property allows us to map the phase coexistence relations onto a 
set on nonlinear equations for the unknown moments of the daughter distribution 
functions. 

We assume that at a certain density po an d composition F (a) the system sep- 
arates into two phases with the densities p\ and arid compositions -^1(0") and 
F2(a). Hereafter the lower index refer to the parent phase and the lower indices 1 
and 2 refer to the daughter phases. At equilibrium these quantities take the values, 
which follows from the phase equilibrium conditions, i.e.: (i) conservation of the 
total volume of the system, (ii) conservation of the total number of the particles of 
each species, (iii) equality of the chemical potentials of particles of the same species 
in the coexisting phases, (iv) equality of the pressure in the coexisting phases. These 
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conditions finally lead to the following set of relations [15|, [l 



F a (a) = F {a) Q a {a; p Q , P i, p 2 ] [F a ]) , (72) 

P 1 (p 1 ;[F 1 ])=P 2 (p 2 ;[F 2 ]), (73) 
F a (a) da = 1, for a = 1 or a = 2, (74) 



where 



n t [F n Po (P2 ~ Pi) {1 - S la + S la exp 

p a Qa\cr;po,PuP2; [F a \) = -. r ( (ex) . , 75 

po- pi - {po - p2) exp {(3Ap^ x >) 

V" = P^ Ps; t^]) - (a, Pl ; [Fx]) , (76) 

PaY exS> is the excess (over the ideal gas) chemical potential of the particle of species 
a in the phase a and [. . .] denote functional dependence. The relation between 
F (cr) and daughter phase distribution function F a (a), i.e., Eq. (1721) . follows from 
the phase equilibrium conditions (i)-(iii). 

Relations (172l) -( J74l) represent a closed set of equations to be solved for the un- 
knowns p a and F Q (a); this set have to be solved for every value of the species 
variable a. However, since thermodynamical properties of the model at hand are 
defined by the finite number of the moments we can map this set of equations onto 
a closed set of 10 + 28iV y + 6N Y {N Y + 1) algebraic equations for p a , C$$ a) , c£" )(a) 
and moments m!^\ m!j?\ rn^ a \ m^°^ a \ rn^f 1 ^ in the two coexisting phases 
(a = 1,2). We have 

poo 

™ ( k ] =Pa rn k (a)F (a) Q a (a, p ; {XJ , {X 2 }) da, k = 0,1, 2, 3, (77) 
Jo 

POO 

ml a) =p a ml a \a)F (a)Q a (a,p ;{X 1 },{X 2 }) da, (78) 
Jo 

POO 

™ { Ki ia) =P« mP} a \a)F (a)Q a (a, Po ;{X 1 },{X 2 }) da, i = 0,1,2, (79) 
Jo 

POO 

^kT ] =pJ m^ {a \a)F (a)Q a (a,p ;{X 1 },{X 2 }) da, j = 0,l, (80) 
Jo 

POO 

m C^X«) =p / m^ a \a)F (a)Q a (a,p ;{X 1 },{X 2 }) da, 1 = 0,1,2, (81) 
Jo 

where K = A,B and {X a } represent unknowns of the problem, i.e. 

{X a } = { Pa ,mi a \ m<?\ mP} a \ m%f a \ m^ )[a) } , a = 1,2. 
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The remaining 2 + 8Ny equations follows from the equality of the pressure in coex- 
isting phases (1731) . 

P 1 (p 1 ;{X 1 })=P 2 (p 2 ;{X 2 }), (82) 

from the set of equations for C^j and C^ 2 °^ (1531) and from the normalizing 
condition (1731) for either phase a = 1 or a = 2, 

/■oo 

/ F (a) Q a (a, p ; {X 1 } , {X 2 }) da = 1. (83) 
Jo 

Solution of the set of equations (|53|) . (1771) - (lB3j) for a given density po and distribution 
function F (a) of the parent phase gives the densities p a and distribution functions 
F a (a) of the two coexisting daughter phases. The coexisting densities at different 
densities of the parent phase po defined binodals, which are terminated when a den- 
sity of one of the phases is equal to the parent phase density po- These termination 
points form the cloud and shadow coexisting curves. These curves intersect at the 
critical point, which is characterized by the critical density p cr = pi = p 2 = po and 
critical temperature T cr . The cloud-shadow curves can be obtained as a special so- 
lution of the general coexisting problem, when the properties of one phase are equal 
to the properties of the parent phase: assuming that the phase a = 2 is the cloud 
phase, i.e. p 2 = po, and following the above scheme we will end up with the same 
set of equations (l53j) . (PT7j) - (l8"3|) . but with p 2 and F 2 (a) substituted by p and F (a), 
respectively. 



IV. RESULTS AND DISCUSSION 



In this section we present our numerical results for the liquid-gas phase diagram 
of polydisperse dipolar hard-sphere fluid at different degree of polydispersity. For 
size distribution function F(a) we have chosen the beta distribution given by 

F{a) = I ^ J — H K -a)H(a- a d ) , (84) 

Y{u)Y{p) [a u -a d ) 



where 



cr-(Td 1-ctq ( 1 n \ o- -a d 
x = , v = — a , p=[- l)u, a = , (85) 



&u - &d Da ' V°"0 / CTu — Vd 

a = (a), D a = ^-1, (a n )= I da n aF(a). (86) 
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Here a u and defined the range of values for a. In our calculations we have chosen 
o u = 1.2947<7o and o"d = 0.85cto- 

In Figs 1-3 the liquid-gas phase diagram for polydisperse dipolar hard-sphere 
fluid at different degree of polydispersity D a = 0.1, 0.2, 0.3 is presented in T* vs p* 
coordinate frame. Here T* = ksT / d 2 ^(ao) and p* = p * erg. We show the cloud and 
shadow curves and critical binodal. In addition for the reference we include Monte- 
Carlo predictions for the critical point and ETPT-CF predictions for the phase 



diagram of monodisperse version of the model [111]. One can see that upon increasing 
of D a the region of the phase instability increases with the critical point shifting to 
the higher temperatures and densities. For the larger values of D a (D a = 0.2, 0.3) 
the low density part of the cloud curve and the high density part of the shadow 
curve become almost vertical. For D a = 0.3 the cloud and shadow curves do not 
intersect. At T* = 0.1715 the cloud curve has a cusp and shadow curve has a jump 
discontinuity. We believe that at this temperature there is a three-phase equilibria, 
when the mother phase is in equilibrium with two phases on two branches of the 
shadow curve, one with slightly lower density and the other with slightly higher 
density, respectively. The cloud and shadow curves for the whole set of values for 
D a are collected in Fig. 4. In Figs 5-7 and Figs 8-10 we present the the average A-size 
La and £>-size L B of the clusters, respectively, along the cloud and shadow curves 
and along binodals for both polydisperse and monodisperse versions of the model. 



For L K we have used expression fl38|) . extended to account for polydispersity, i.e. 



J K 



4- + 

daF(a 



[i + 4»] [i + k|> 



4 -(!-«*») 



dvF(a) =- ■ - - -= i - ' J } . (87) 



[l + 4»] + 

From these figures one can see that clusters of the larger sizes occurs in the liquid 
phase, in comparison with the gas phase. Decrease of the temperature causes in- 
crease of the cluster sizes in the liquid phase and decrease in the gas phase. With 
the increase of polydispersity along the cloud curve do not change much. How- 
ever corresponding changes along the shadow curve is more substantial, here the 
cluster sizes increases with polydispersity increase. In Figs 11-13 we show the ratio 
of A— and 5-sizes La/Lb- As one would expect La is substantially larger then 
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Lb, thus prevailing is a chain structure of the clusters formed, with chains mutu- 
ally connected via 5-bonds. This difference in La and Lb is larger in the liquid 
phase and increases with polydispersity increase. Similar, as before, the temper- 
ature decrease causes La/ Lb decrease in the gas phase and increase in the liquid 
phase. In Figs 14-16 and Figs 17 and 18 distribution functions of the shadow curve 
and on the critical binodal at different temperatures are demonstrated. According 
to these figures the larger size particles always fractionate to the liquid phase and 
smaller particles prefer to move to the gas phase. With increase of D a and decrease 
of the temperature this fractionation effects becomes more pronounced. Finally in 
Fig 19 we show distribution functions of the two branches of the shadow curve at 
T* = 0.1715 and mother phase distribution function. 



V. CONCLUSIONS 



In this paper we propose extension of our TPT-CF approach to account for 
several associating sites with the possibility of each site to be multiply bonded. 
The theory is applied to study the liquid-gas phase behavior of polydisperse dipolar 
hard-sphere fluid. We present a full phase diagram, which include cloud and shadow 
curves, binodals and distribution functions of the coexisting phases. 



VI. APPENDIX A 

Orientationally averaged Mayer function ]kk{j'- i <J\ (J 'i) were fitted empirically as 
a sum of Yukawa-like terms 

N -z (m) (T)(r-cr- ■) 

fm(r, aw) = 4"° T ) 4"° fo> T ) 6 " M frr . 13 . (88) 

n=l Z n (T) T 

where f\ = Jaa and f 2 = Jbb, hj denotes the species of the particles, cr^ = 
(<7j + <7j) /2, N = 6, and A n (cr, T) are parameters that depend on particle size 
and temperature and were fitted by "polynimally-exponential" functions of different 
forms for first and second integrals. 

The fitting was performed for following range of parameters: a G [0.85, 1.2947] 
and T e [0.13,0.2]. The functional dependence of A n (<7j,T) and z n (T) was chosen 
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differently for first and second integrals. For the first integral 

4 1} T) = a® (T) + a% (T) x + a% (T) x 2 + a£> (T) x 3 



+ a^T) e a S (T) * +a »^^ (89) 

where x = (cr — o" m j n ) /(Acr), <7 min = 0.85, Act = 1.2947 — cr m j n = 0.4447 and 
a n\ {T) ,—, Un,9 (T) are given below 

<! (T) = A* (&Sl« + 6 SU + 6 £U 2 ) > ^ or * = 2 > 3 ' 4 

<! CO = (&£i,i + O + &£U 2 + O 3 ) . /or 1 = 5,6,7,8,9, (90) 

where y = and T min = 0.13. 

The temperature dependence of zjP (T) reads 

4 1} CO = + + 4V + 4V) • (91) 

(2*1 f2) 

For the second integral the functional dependence of A n (cr, T) and % (T) is 
following: 

4 2) (T) = (a$ + + 4V + 4^^) (92) 

where 

aS = A" (&S, 2 + O + 6£lieC»), /or i= 1,2,3,4 (93) 

and 

41 = (d + O + &S,3l/ 2 + &S,4^' 52/ ) , /or < = 5, 6, 7, 8, 9. (94) 

The x and y are the same as for the first integral. 

The fitting procedure consisted of finding suitable b n ^j and u n j parameters for 
irst and second integrals by means of differential evolution optimization algorithm 
I?]) Q|. The objective function in both cases was a sum of square deviations of the 
area under f m (r,(Ji(Tj) as a function of r, and contact value f m (r = , (JiOj) , from 
theirs fitting representations for different values of <7j, aj and T. The deviations in 
the objective function were calculated for ten different temperature values (uniformly 
distributed from 0.13 to 0.2) and twenty different a values (also uniformly distributed 
from 0.85 to 1.2947). 
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VII. APPENDIX B 
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FIG. 1: Predictions of the ETPT-CF for the phase diagram of polydisperse dipolar 
hard-sphere mixture including cloud and shadow curves (as labeled), and critical bin- 
odal (dashed line) at D a = 0.1 in p* vs T* coordinate frame. Dotted line and empty 
square represent ETEP-CF binodal and MC [? ] critical point of monodisperse dipolar 
hard-sphere fluid, respectively 
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FIG. 2: The same as in Fig.l at D a = 0.2. 




i i i i i r 

cloud shadow 




24 



0.19 



1 1 1 ~7\ I 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 
P 



FIG. 3: The same as in Fig.l at D a = 0.3. 
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FIG. 4: ETPT-CF cloud and shadow curves for polydisperse dipolar hard-sphere mixture 
at D a = 0.1 (solid lines), D a = 0.2 (dashed lines) and D a = 0.3 (dotted lines). 
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FIG. 5: Average yl-size of the clusters along the cloud and shadow curves (as labeled), 
along the critical binodal (dashed line) for polydisperse dipolar hard-sphere fluid at D a = 
0.1 and along the binodal for monodisperse version of the model (dotted lines). 
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FIG. 7: The same as in Fig.8 at D a = 0.3 
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FIG. 8: Average S-size of the clusters along the cloud and shadow curves (as labeled), 
along the critical binodal (dashed line) for polydisperse dipolar hard-sphere fluid at D a = 
0.1 and along the binodal for monodisperse version of the model (dotted lines). 
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FIG. 9: The same as in Fig.ll at D a = 0.2. 
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FIG. 10: The same as in Fig.ll at D a = 0.3 
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FIG. 11: The ratio of the average A and B sizes of the clusters along the cloud and shadow 
curves (as labeled), along the critical binodal (dashed line) for polydisperse dipolar hard- 
sphere fluid at D a = 0.1 
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FIG. 12: The same as in Fig. 14 at D a = 0.2 
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FIG. 13: The same as in Fig. 14 at D a = 0.3 
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FIG. 14: Distribution functions of the gas (dotted lines) and liquid (dashed lines) phases 
along the shadow curve for T* = 0.16, 0.15, 0.14 and mother phase distribution function 
(solid line) at D a = 0.1. With the temperature decrease distribution functions of the 
liquid phase shifts in the direction of larger a and distribution functions of the gas phase 
shifts in the direction of smaller a. 
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FIG. 15: The same as in Fig.17 for T* = 0.17, 0.16, 0.15, 0.14 and D a = 0.2 
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FIG. 16: The same as in Fig.17 for T* = 0.18, 0.17, 0.16, 0.15, 0.14 and D a = 0.3 
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FIG. 17: Distribution functions of the gas (dotted lines) and liquid (dashed lines) phases 
along the critical binodal for T* = 0.16, 0.15, 0.14 and mother phase distribution function 
(solid line) at D a = 0.1. With the temperature decrease distribution functions of the liquid 
phase slightly shifts in the direction of larger a (almost coinsiding with the mother phase 
distribution function) and distribution functions of the gas phase shifts in the direction of 
smaller a. 




FIG. 18: The same as in Fig.20 at D a = 0.2 
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FIG. 19: Distribution functions on a shadow curves at T* = 0.1715 for three phases at 
equilibria for D a = 0.3 
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